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Abstract 

Electromagnetic waves, solving the full set of Maxwell equations 
in vacuum, are numerically computed. These waves occupy a fixed 
bounded region of the three dimensional space, topologically equiva- 
lent to a toroid. Thus, their fluid dynamics analogs are vortex rings. 
An analysis of the shape of the sections of the rings, depending on 
the angular speed of rotation and the major diameter, is carried out. 
Successively, spherical electromagnetic vortex rings of Hill's type are 
taken into consideration. For some interesting peculiar configurations, 
explicit numerical solutions are exhibited. 
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1 Introduction 



We start by recalling that the classical set of Maxwell equations in vacuum 
has the form: 

^ = c^curlB (1) 

dt ^ ' 

divE = (2) 
^ = -curlE (3) 

divB = (4) 

where c is the speed of light and the two fields E and cB have the same 
dimensions. 

It is customary to introduce the potentials A and $ such that: 

B = - curlA E = - - V$ (5) 

c c dt ^ ' 

By this assumption, equations ([3]) and (jH) are automatically satisfied. In 
addition, we require that the vector and the scalar potentials are related by 
the following Lorenz gauge condition: 

1 5$ 

divA + = (6) 

c ot 

Furthermore, one can deduce the two wave equations: 



- c^AA = (7) 



^2 



A$ = (8) 



We are concerned with studying, from the numerical viewpoint, the de- 
velopment of a solitary electromagnetic wave, trapped in a bounded region 
of space having a toroid shape. In the cylindrical case (equivalent to a 2- 
D problem), and partly for 3-D problems, this analysis was proposed and 
carried out in Ref. [Ij. There, the goal was to simulate stable elementary 
subatomic particles by means of rotating photons. Here, we would like to 
continue the discussion of the 3-D case, because the subject might be of 
more general interest. Since exact solutions are only available in special sit- 
uations, we shall make use of finite element techniques to approximate the 
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model equations. In order to determine solutions of the vector wave equa- 
tion, confined in appropriate steady domains, we will perform an in-depth 
analysis of the lower spectrum of a suitable elliptic operator, in dependance 
of the shape and magnitude of the regions. In particular, we will be con- 
cerned with those domains realizing the coincidence of the fourth and the 
fifth eigenvalues of the differential operator. The reasons for this choice will 
become clear to the reader as we proceed with the investigation. 

The structures we consider in this paper display a strong analogy with 
fluid dynamics vortex rings (see for instance Ref. [2J, [3j, [4J). For this reason, 
the techniques we apply here may be useful to get additional results in the 
study of the development of fluid vortices. As a matter of fact, peculiar 
configurations are going to be presented and discussed, opening the path to 
an interesting scenario for future extensions. 



Let us quickly review the results obtained in Ref. jl], concerning waves 
rotating around the axis of a cylinder. The magnetic field is oriented in the 
same direction as the z-axis, so that the electric field lays on the orthogonal 
plane. The variables are expressed in cylindrical coordinates (r, z,^). The 
solutions, however, will not depend on z. 

Several options are examined in Ref. [Ij, here we show the most signifi- 
cant one, that will be used later to construct the toroid case. We denote by 
uj a positive parameter which characterizes the frequency of rotation of the 
wave: v — cuj/ti. Up to multiplicative constant, the two vector fields are 
given by: 



for < < 27r, < r < 5o/cj and any z. The fields in (|9j) are defined on a 
disk of radius ^q/cj, whose size is inversely proportional to the angular speed 
of rotation. In Ref. [1] , figure 5.7, the reader can see the displacement of 
the electric field for t = (some on-line animations can be viewed in Ref. [5j, 
rotating photons). 

In (|9j) we find the Bessel function Jj^ (see, e.g., Ref. [6j) that solves the 



2 The Cylindrical Case 





(9) 



3 



following differential equation: 



J'ki^) + ^ - + M^) - (10) 

It is also useful to recall that Bessel functions are connected by the relations: 

+ = Jk-iix) (11) 

J,+i(x) = - J,_^ix) (12) 

The quantity Sq ^ 5.135622 turns out to be the first zero of J2. In this 
way, for r = Sq/oj, the components Ei and are zero. Note also that E 
and B vanish for r = 0, since Jk{x) decays as for x ^ 0. The choice 
A: = 1 is not permitted because it does not allow to prolong with continuity 
the fields up to r = 0, although, one could take into consideration solutions 
for k integer greater than 2: 

E = ^Qg(^^^^ _ ^^^^ J'k{(^r) sin(ccjt - kcj))^ 

B = -fo, Jk{cjr)cos{ccjt-k(f)), o) (13) 



c 

The idea is to simulate a fc-body rotating system in equilibrium. This some- 
how explains why the case A: = 1 is not going to produce meaningful solu- 
tions. 

For k — 2^ the electromagnetic fields in ([9]) are generated by the following 
potentials: 

A = — — (^J-^iur) siiiicut — 20), , Jz{uor) cos(ccjt — 20)^ 
$ = - — J2(cjr) cos(ccjt - 20) 

UJ 

satisfying the Lorenz condition ([6]) and the equations ([7j)-([8]). The reader 
can check this by direct differentiation. We observe that Ai and have a 
phase difference of 45 degrees, since: sin(ca;t — 20) = cos(ccL;t — 2(0 + 7r/4)). 
Denoting by a(r) = J^{ujr) the term depending only on the radial variable, 
we easily discover that (see (|TOl) for /c = 3): 

1^__9_\ _ 

dr^ r dr r'^ I 
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We can introduce a new function w = -^a -\- ^a. A straightforward compu- 
tation, using relation (fT4l) and its derivative, brings to: 

/ d2 Id 4 \ _ / /// a^' 10a' 18a 

\ dr^ r dr r'^ J \ r 



3 / a' 9a\ 9 / / 3a 

- a H 

r V r r 



+ = c.2(a' + :^) = (15) 



By scaling the interval in order to impose the boundary conditions 
w{0) = and w{l) = 0, the eigenvalue problem (IT5]l takes the form: 

- -[d;;:2 + -^,--2)-' - ^0- (16) 

where L is a positive-definite differential operator. Such an operator is 
the same as the one we would obtain from the Laplacian, after separation of 
variables in polar coordinates, by imposing homogeneous Dirichlet boundary 
conditions on a disk (on the plane (r, 0)) of radius equal to 1. In this 
circumstance, the eigenvalue: 

X = 6^ ^ 26.37461 (17) 

has double multiplicity. The corresponding eigenfunctions are orthogonal in 
L^(fi) and show a phase difference of 45 degrees (see figures 1 and 2). 



A 


multiplicity 


5.78318 


1 


14.68197 


2 


26.37461 


2 


30.47126 


1 


40.70646 


2 


49.21845 


2 


57.58294 


2 


70.84999 


2 



Table 1: Eigenvalues counted with their multiplicity of the Laplacian on a disk of 
radius 1, with homogeneous Dirichlet boundary conditions. 

In other words, by computing the spectrum of the Laplacian on a disk of 
radius 1 (see table 1), the one given in (fT7|) is the common eigenvalue of the 
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fourth and the fifth eigenfunctions. Note that there are no eigenvalues with 
multiphcity greater than 2. We also note that $ = at the boundary of Q. 
These simple observations will be of primary importance in the discussion 
to follow. 



a 3 1 ^0^*t 

Figure 1: Signature of the first 8 eigenfunctions of the Laplace's equation on a 
disk. The first eigenfunction (wi) does not change sign. The two successive ones 
{w2 and ws) display a phase difference of 90 degrees. Then, we have and 
with a phase difference of 45 degrees. The next one (wq) is a single multiplicity 
eigenfunction, while the last ones {wr and ws) have a phase difference of 30 degrees. 



Figure 2: Solutions of the wave equation on a disk are obtained by linearly combin- 
ing independent eigenfunction having the same eigenvalue. For example, a complete 
rotation is simulated through the sequence: W4,W5, —W5,W4,W5, —w^. 
Intermediate situations are obtained by taking W4 cos cty/X + sin cty/X, with 
< cty/X < 27r, where A is the common eigenvalue. 

3 The Toroid Case 

We are ready to study the 3-D case. In cylindrical coordinates (r, 0), we 
consider the potentials: 

A = 1 [a{t, r, z), bit, r,z),0) ^ = - + 7 + (18) 



6 



where a and h are functions to be computed, while A and B are their primi- 
tives with respect to the time variable. In this context, the Lorenz condition 
([6j) turns out to be automatically satisfied. The functions a and h will de- 
scribe the time evolution, in a certain region of the plane (r, z), of a 
rotating wave similar to the one studied in the previous section (differently 
from the previous case, is now vertically oriented). Since there is no de- 
pendance on 0, the solution is automatically extended to a toroid S, having 
section fi, with the axis parallel to the z-axis. 

From the potentials we deduce the electromagnetic fields (see ([5])): 

^ r I da I dh \ ^ I r ^ dh da\ 

\ dr dt^ dz dt^ J V ' ' dr dz J 

We now require E and B to satisfy the whole set of Maxwell equations (in 
alternative, we get the same conclusions by imposing ([7]) and (jS])). The 
result is the following system: 



1 d'^a d'^a d f da 



a 



+ ^Air + z) (20) 



dt'^ dz'^ dr \dr r 



1 d^h _ d'b d^b Idb 
dt'^ dz'^ dr'^ r dr 



We can couple a and b through the boundary condition: 

da db 
dz dr 



(22) 



which amounts to ask B = c ^ curl A = at the contour of ft. The 
field E is tangential to the same boundary. As a matter of fact, thanks to 
([ISl)-(Eni)-(EI]), one can write: 

E - -l^(r— -r—] o] (23) 

\dz\dz dr r dr\ dz dr J 

which means that E is orthogonal to the gradient of r(^-^A — ■§pB^. Since 

(f22l) is valid for any t, such a gradient is orthogonal to the boundary of fi, 
showing that E is tangential. 

Note that now the set Q is not going to be a circle. As done in Ref. [1], 
section 5.4, we set y — r — rj. The domain will be centered at the point 
(77, 0), where 77 > is large enough to avoid intersection of with the z-axis. 
The quantity 277 is related to the major diameter of the toroid. 
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By differentiating the first equation with respect to z and the second one 
with respect to one gets: 



1 d^u d^u ^ d fdu ^ u 
dt'^ dz^ dy \dy y + rj 



1 d^v d^v d f dv V 
. dt'^ dz'^ dy \dy y + r] 



(24) 



where we took u — and v — 
thermore, with the substitutions u — u^y + 77, v 



^ 6, with u — V 



dy 



on dVt. Fur- 
v^y + 77, one obtains: 



( 1 d'^ii 


_ 92n 


5y2 


3 n 




dz^ 


4 (y + r?)2 


1 




9y2 


3 V 


I C2 5*2 




4 (y + 17)2 



(25) 



to be solved in fi, with the boundary condition: u — v — 0. 

Similarly to the case examined in section 2, we look for solutions associ- 
ated to the second mode (A: = 2 in (1T3]1 ) with respect to the angle of rotation. 
According to the captions of figures 1 and 2, the functions u and v should 
have a phase difference of 45 degrees (this property is now qualitative, since 
is not a perfect disk). We fix the area of to be equal to tt. For 77 tending 
to infinity, the set converges to a circle of radius 1 and the electromagnetic 
fields coincide with those given in (|9j) for = ^q. 

We proceed by introducing a new unknown w — u — v and by taking the 
difference of the two equations in (125]) . Then, we can get rid of the time 
variable and pass to the stationary eigenvalue problem: 



Lw 



w 



dz^ 
on 



dy'^ 



+ 



w 



4 + 7/)2 



Xw in Q 



(26) 



Thus, we would like A > to be an eigenvalue of the positive-definite differ- 

+ |(y + 77)"^, with homogeneous Dirichlet 



0^2 



ential operator: L 

boundary conditions. In addition, since we want two independent eigen- 
functions with a difference of phase of 45 degrees, the multiplicity of A must 
be equal to 2. This implicitly defines More precisely, we require that 
A = A4 = A5, where A4 and A5 are the eigenvalues of the fourth and the fifth 
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eigenfunctions, W4 and w^, of L on ft. In order to preserve energy, these 
eigenfunctions will be normalized in L^(fi). In this way, by setting: 

u{t^y^z) — v{t^y^z) — W4^{y^z)^\n[ct\f\) + w^{y ^ z) cos{ct\f\) (27) 

one gets a full solution of the system (l25l) . We can finally recover the elec- 
tromagnetic fields by the expressions (see (l23ll ): 



E 



-{-{'- 

a/a V dz \ 



-W4 cos( + sm( 



-1 d 



VyTvdy 



{-W4 cos C + '^5 sin C)\/y + r/ ) , 



with ( = ct\f\. As in figure 2, by varying the parameter t, we can simulate 
a rotating object. 

Of course, we could also take into consideration the case k > 2 (see the 
2-D analog (|T3l) ) by studying the behavior of eigenvalues with higher mag- 
nitude. The search, in this case, should be addressed to the determination 
of couples of independent eigenfunctions sharing the same eigenvalue. We 
think this is a viable option, although, in order to maintain the discussion 
at a simple level, we will not analyze this extension. 

We now define H = x [0, 27r[, in order to get a 3-D solution not de- 
pending on (j). In fluid dynamics this structure is known as vortex ring (see 
for instance Ref. [2J, section 7.2, or [3]). On the surface of H, field E is 
tangential and B is zero. At every point inside S, the time average dur- 
ing a period of oscillation, of both E and B, is zero (see the animations in 
Ref. [5]). 

We would like to know more about the shape of S. Considering that not 
all the sets are such that A4 and A5 are equal, we can use this property in 
order to determine the right domain. This problem admits however infinite 
solutions. In the next section we discuss how to find numerically some 
suitable configurations. 



4 On the Optimal Shape of Q 

We first observe that, among the sets with fixed area equal to tt, the circle 
of radius 1, minimizes the eigenvalues of the Laplacian with homogeneous 



9 



Dirichlet boundary conditions (see for instance Ref. [7J). Due to symmetry 
arguments, all the eigenvalues related to the angular modes have multi- 
plicity 2 (see table 1 and figure 1). The eigenfunctions are supposed to 
be orthogonal and normalized in L^(fi). Then, the normal derivative of 
the eigenfunction related to the lowest eigenvalue, has constant value on 
the boundary of the disk. Similarly, the sum of the squares of the normal 
derivatives of the eigenfunctions related to the second and the third eigen- 
values, is constant along the boundary. The same is true for all the couples 
of eigenfunctions relative to other eigenvalues with double multiplicity (the 
fourth and the fifth, for example). 

We now replace the Laplacian by the new elliptic operator L — — 

^ + |(y+77)~^ (see (l26]l ). with homogeneous Dirichlet boundary conditions. 
We are concerned with finding a set Q, with area equal to tt, such that the 
fourth and the fifth eigenvalues, A4 and A5, of L are coincident. Thanks to 
(|27|) . this allows us to determine solutions of a wave- type equation, rotating 
inside with an angular velocity proportional to a/A, where A = A4 = A5. In 
table 2, we report the eigenvalues of L on the disk centered in (y, z) = (0, 0) 
and radius equal to 1. If < 77 < 1 the differential operator is singular inside 
(the corresponding toroid region S has no central hole). Therefore, we 
take 77 > 1. The case 77 = 1, where the axis y — —rj (or, equivalently, r = 0) 
touches at the boundary, is still admissible. The table shows the results 
for different values of 77. For 77 large, the fourth and the fifth eigenvalues are 
quite similar, due to the fact that the term |(^ + 77)~^ becomes small. For 
77 close to 1, the situation is not too bad, anyway. This is especially true for 
the highest eigenvalues, since the corresponding eigenfunctions, due to the 
boundary constraints, decay faster near the border of Vt. Therefore, they do 
not "feel" too much the presence of the term |(^ + 77)"^. 







rj = 1.2 


ri = 1.5 


17 = 2 


Ai 


6.90 


6.46 


6.18 


5.99 


A2 


15.71 


15.33 


15.08 


14.90 


A3 


16.75 


15.68 


15.19 


14.93 


A4 


28.06 


27.27 


26.86 


26.63 


A5 


28.36 


27.35 


26.88 


26.64 



Table 2: The first 5 eigenvalues of the operator L on a disk of radius 1, for different 
values of rj. 
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For any fixed 77, we would like to adjust the shape of ft with the help 
of some iterative procedure. The idea is to correct the boundary in order 
to fulfill a certain condition at the limit. Before getting some interesting 
answers (see later), we tried unsuccessfully different approaches. Let us 
discuss the one that looked more promising. We consider the substitution: 



w{y,z) 



In this way f[26|) becomes: 



w{s, z) 



where s — {y ^ riY 



(29) 



(^4 ^ 
w — on 50 



\w in O 



(30) 



where now the operator is a Laplacian with a variable coefficient. Concerning 
the operator = ~ for /3 > constant, the domain O of fixed 

area equal to tt, optimizing the set of eigenvalues, is an ellipse. As in the case 
of the circle (/3 = 1), we have infinite eigenvalues with double multiplicity. 
In addition, if w is the first eigenfunction of L^, the following relation is 
easily checked: 



constant on 90 



(31) 



The idea is to generalize the above relation in the case of variable coefficients. 
Hence, let us suppose that w is an eigenfunction corresponding to the first 
eigenvalue in ([501) . we may require that: 



+ 



/dwy 



constant on dQ 



(32) 



For points where s is large, we expect the curvature of the boundary of O 
to be high, and viceversa. The resulting H is a kind of doughnut, a bit 
flattened on the internal side. Due to ([29]), in terms of w one has: 



(y + v) 



{y + v) 



w 



+ 



dw ^ ^ _^ ^ dw ^^ 2 



2y/y^r] dyJ \dz J 



\dyJ 

where we used that it; = on 50 



\dzJ 



constant 



(33) 
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Relation f[33|) is the "target" we would like to reach by implementing 
our iterative method. To this end, starting from the circle of radius 1, 
we compute the eigenfunction (normalized in L^(fi)) corresponding to the 
lowest eigenvalue in ([261) . We then compute its normal derivatives on dQ 
and correct in order to enforce condition (f551) . For points belonging to 
the boundary, the updating procedure is performed as follows: 

Pnew = Fold + ^^^N (34) 

where N is the outer normal derivative, is a relaxation parameter and £ 
is the difference between the normal derivative evaluated in Fold the 
average of the normal derivatives computed on the boundary of the current 
fi. After every correction, the area of is set to tt, with the help of a linear 
transformation. At the limit, we would like to have = 0. 

In order to compute eigenvalues and eigenfunctions we used a finite ele- 
ment code. In particular, we implemented MODULEF (see Ref. [5] and [H]) 
with P2 elements. The grid was fine enough to have reliable results up to the 
second decimal digit (at least for the first eigenvalue). The solver is based 
on a QR algorithm. The points to be updated through the iterative method 
(l34]) are the vertices of the triangles belonging to dQ. The normal derivative 
at these points is the average of the normal derivatives at the mid-points 
of the two sides of the contiguous triangles. After the entire boundary has 
been modified, a brand-new grid is generated in fi. The performances of 
this technique are far from being optimal, but their improvement is not in 
the scopes of the present paper. 

Unfortunately, despite all the efforts, the method did not want to con- 
verge. A possible explanation is that relation (f33|) cannot be realized, be- 
cause it is wrong from the theoretical viewpoint. The iterative technique 
(l34]) was however a good starting point to try successive variants, and, fi- 
nally, we had success with the following scheme: 

Fnew = Fold + ^|A4-A5| (n^^^ {v old + vh^"^) ondn (35) 

where N = {rSy\n^^^) is again the outer normal derivative. This time, 
the fourth and the fifth eigenvalues of L explicitly appear in the correcting 
term. Inspired by (133]) . the updating on the boundary is not uniform in the 
variable but the second component is weighted by the function y + 77. 

After suitably adjusting 0, convergence is obtained in a few iterations. 
The result is a domain yielding A4 = A5. However, as anticipated, there 
are infinite other domains with this property. The most interesting ones 
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should be those having a distribution of the eigenvalues as lower as possible, 
and making the shape of as rounded as possible. We do not know if 
what we got by iterating ([35l) can be considered optimal with this respect. 
Nevertheless, it is important to have shown that the condition A4 = A5 can 
be actually achieved (at least numerically). In table 3, one finds the modified 
eigenvalues for various 77 (compare with table 2). Note instead that, for 77 
close to 1, A2 and A3 are not in good agreement. 

For 77 = 1 and 77 = 1.2, the corresponding domains are displayed in 
figure 3. Even for 77 close to 1, they are not too far from a circle. In figure 
4, for 77 = 1, we provide the plot of the fourth and the fifth eigenfunctions, 
actually showing a phase difference of 45 degrees (see also figures 1 and 2). 
For 77 > 2, is practically a circle. 

The analysis carried out in this section may help to understand the 
structure of fluid dynamics vortex rings. These are incredibly stable con- 
figurations, that have been widely studied, under different aspects. Among 
the vast literature we just mention for instance the papers: [10], [11], [4], 
fl2] . In addition, in Ref. |13j, "optimal" sections of vortex rings are dis- 
cussed. However, in the fluid dynamics case, the situations where the major 
diameter is relatively small, bring to more pronounced deformations of the 
sections. 





T] = 1 


7] = 1.2 


T] = 1.5 


r? = 2 


Ai 


6.85 


6.44 


6.17 


5.98 


A2 


15.76 


15.32 


15.04 


14.89 


A3 


16.59 


15.64 


15.18 


14.90 


A4 


28.08 


27.24 


26.83 


26.60 


A5 


28.08 


27.24 


26.83 


26.60 



Table 3: The first 5 eigenvalues of the operator L on the modified domain, for 
different values of r]. 

We are not able to discuss the stability of our electromagnetic vortices. 
Theoretically, using the newfound 3-D solutions, one should build the cor- 
responding electromagnetic stress tensor and put it on the right-hand side 
of Einstein's equation, as illustrated in Ref. \T\. Then, one has to solve this 
nonlinear system in order to determine the metric tensor. Finally, one should 
check that the quasi-circular orbits (fi is not a perfect circle) followed by the 
light rays are actually geodesies of such a space-time environment. Note that 
this study depends on the parameter cj, connected to the speed of rotation 
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of the wave. For uj large, the size of is small and the frequency is high. 
We expect one or more situations of equilibrium, where the gravitational 
setting and the centrifugal effect of the spinning solitons are compensated. 
This would mean that only particular geometries are admitted, specifying 
exactly the shape and the size of the toroid regions. 




Figure 3: Toroid sections (solid line) compared to the circle of radius 1 (dashed- 
dotted line), for 77 = 1 and 77 = 1.2. The vertical line is the toroid axis of symmetry. 




Figure 4: The forth and the fifth normalized eigenfunctions corresponding to the 
case 77 = 1. Even if is not a perfect circle, they actually display a phase difference 
of 45 degrees. Each eigenfunction has two positive (white) and two negative (black) 
bumps. 
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Of course, what we just said above turns out to be very hard to prove, 
both theoreticahy and computationahy. Note that we are also omitting a 
stationary component of the electromagnetic fields, introduced in Ref. [Ij, 
which should provide charge and mass to the particle model. We do not 
insist further on this topic and we address the reader to Ref. pLj for additional 
information. 

5 HilPs Type Vortices 

The next step is to examine what happens outside the spinning toroid. If 
we were dealing with a fluid vortex, due to viscosity, we might expect the 
formation of other external vortices developing at lower frequency. This 
circumstance can be observed in tornados or typhoons (see for instance 
Ref. [14J). In the theory presented in Ref. [Ij, an electromagnetic vortex, 
through a mechanism still to be clarified, captures the surrounding electro- 
magnetic signals and generates a series of encapsulated shells vibrating with 
decreasing frequencies. These shells could be responsible for the quantum 
properties of matter. 




Figure 5: Successive layers of encapsulated vortex rings. In the first picture, the 
toroid regions have a common central hole. In the second one, the central hole is 
reduced to a vertical segment (Hill's vortex). Of course, intermediate situations 
combining some shells of the first type, successively embedded in spherical vortices, 
could be taken into consideration. 
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Basically, there are two possible ways in which external ring vortices 
may develop. As the first picture of figure 5 shows, we may have a series 
of successive toroid structures, where all the fluid stream-lines pass through 
a common central hole. This situation is difficult to analyze, since we have 
no idea of the shape of these regions and the location of the primary vortex 
inside them. In practice, there are too many degrees of freedom to work 
with. 

The other situation (second picture of figure 5) is more affordable. It 
represents an Hill's spherical vortex (see Ref. [2J, section 7.2, and Ref. fT5] 
for some computational results), successively surrounded by other spherical 
layers. Thus, we now know exactly the shape of these structures. We just 
have to find the location (and the relative size) of the primary toroid vortex 
inside the most internal sphere. This research will lead us to interesting 
conclusions. 




Figure 6: Sections of spherical vortices. 



We start by considering the domain A in figure 6. We set the radius equal 
to 1, so that the area is 7r/2. We solve iuQ = A the eigenvalue problem (l26]l . 
involving the operator L = ~ ^ + {r — y -\-ri with 77 = 0). In the 

first column of table 4, we report the approximated eigenvalues, obtained 
by discretization with the finite element method. As the reader can notice, 
the fourth and the fifth eigenvalues are not coincident. Therefore, for a 
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perfect spherical vortex, we have no chances to obtain solutions of the time- 
dependent wave equation. In alternative, we could modify a bit the shape 
of the domain A or accept solutions whose stream-lines are not stationary. 
This is not however the path we would like to follow. If the rotation inside 
Q is generated by an internal spinning toroid, we may cut out a hole (as in 
the domains B and C) and see what happens to the spectrum of L. For 
simplicity, the hole will be a circle of a certain radius, suitably placed in a 
specific spot (although we know from section 4 that such a hole is slightly 
deformed). Playing with the size and the location of the circle, we look for 
situations in which the fourth and the fifth eigenvalues (A4 and A5) of L are 
the same. 

Surprisingly, this research seems to have a finite number of solutions. 
According to table 4 (columns 2 and 3), we have coincidence of A4 and A5 in 
two particular cases. In the first one, the center of the small circle is placed 
at point (^,0) and the radius is approximately equal to tb = .0505. In the 
second one, the center is at (|,0) and the radius is approximately equal to 
rc — .00520. With the exception of these cases, by varying the magnitude 
and the position of the internal circles, we always found A4 7^ A5. In our 
experiments we did not try anyway all the possible configurations, thus we 
do not exclude the existence of other significant settings. 





domain A 


domain B 


domain C 


domain D 


Ai 


20.24 


30.23 


21.64 


0.63 


A2 


33.29 


34.16 


33.32 


0.99 


A3 


48.94 


54.82 


49.09 


1.45 


A4 


59.82 


68.15 


67.17 


1.99 


A5 


67.11 


68.15 


67.17 


1.99 



Table 4: The first 5 eigenvalues of the operator L, based on the four domains given 
in figure 6, with homogeneous Dirichlet boundary conditions. 

We provide in figure 7 the discretization grids. We show in figures 8 and 
9, the time evolution of the rotating waves, obtained from expression ([271) 
for equispaced values of t with c\/~\t G [0,7r] (half-cycle). In the first case, 
the spherical flow is chained to the inner toroid. In the second case, the flow 
avoids the toroid, still remaining inside the spherical region. The sequence 
looks rather complicated. With the help of some imagination, one can see 
two ant i- clockwise rotating waves, circulating independently in the lower 
and the upper quarters. They have different phases, so that the positive 
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(white) bumps and the negative (black) ones, alternately merge to form a 
single protuberance situated near the center. The pictures only show half 
of the cycle. Then, the sequence restarts with the two colors interchanged. 
The corresponding animations can be found in Ref. [5j (click related papers). 

We point out once again that we are dealing with electromagnetic waves. 
Therefore, we should ask ourselves what happens to the vector fields. Our 
plots actually show the evolution of the function '§^b — -^a, where A = 
c~-'^(a,6, 0) is the vector potential. Using (|T9]1 . (EHj) and (l28]l . one then 
computes E and B. We discover that, at each point, when time passes, the 
tips of the arrows of the electric field turn around, describing elliptic orbits 
(see the animations in Ref. [5j in the case of a circle). The average of E 
during a cycle is zero. The frequency of rotation is globally the same, but, 
depending on the point, it is associated with a different phase, so that the 
general framework seems quite unorganized. 




Figure 7: Discretization grids for the domains B and C. The mesh is finer near 
the vertical axis since the differential operator L is singular there. The degrees 
of freedom (internal and boundary nodes) are 2324 and 1968, respectively. In the 
right-hand side picture, the hole is very small. 



18 




Figure 8: Rotating wave in a Hill's spherical vortex: case of the domain B. The 
sequence is referred to half cycle and terminates with the colors inverted. The 
small inner circle is spinning at a frequency more than 12 times greater, dragging 
the spherical wave along circulating paths. 
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Figure 9: Rotating wave in a Hill's spherical vortex: case of the domain C. The 
tiny inner circle (hardly visible in the pictures) is spinning at a frequency about 
120 times greater. The behavior is qualitatively the same as in figure 8, but now 
the wave do not pass between the circle and the vertical axis. 
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From the size of the hole, we can estimate the frequency of rotation of 
the internal toroid regions. Such a frequency is proportional to the square 
root of a suitable eigenvalue. We take the value \/26.37 ^ 5.135, that, ac- 
cording to table 1, is the square root of the eigenvalue corresponding to the 
fourth and the fifth eigenfunctions, for a circle of radius 1. Afterwards, by 
scaling, we get the frequencies on the reduced circles: 5.135/r^ ^ 101.68 
and 5.135/rc ^ 987.5. Thus, the ratio between the frequencies of the in- 
ner spinning rings and the spherical vortices are given by (see table 4, the 
last line of the second and the third columns): 101.68/a/68.15 ^ 12.31 and 
987.5/\/67.17 ^ 120.48, respectively. Note that these numbers may be af- 
fected by rounding errors (in particular, we expect an error bound of ±0.005 
in the computation of the eigenvalues), therefore they should be taken with 
a little caution. Of course, more trustable results can be obtained with a 
finer mesh. 

Finally, considering the second picture of figure 5, we can study the next 
external spherical shell, whose section is given by the set D of figure 6. Let 
us suppose that the radius of the inner circumference is equal to 1, and let us 
find the outer radius r/^) in order to have A4 = A5. From the experiments we 
deduce ro ^ 5.80. This gives A4 = A5 1.99 (see the last column of table 
4). The number of degrees of freedom, including internal and boundary 
points, is 1857. The two corresponding eigenfunctions are shown in figure 
10. Their time evolution is very similar to the one of figure 9. 

We can confirm and improve these last results by explicitly determining, 
in terms of classical orthogonal basis, the solutions of the wave equation. 
Due to the simplicity of the domain we can actually use separation of 
variables in spherical coordinates (r, 0, 0). Note that here the variable r has 
a different meaning: before it was the distance from the axis (0, z, 0), now it 
is the distance from the center (0,0,0). According to Ref. [6| and Ref. [1], 
p. 15, we can get solutions to the spherical vector wave equation by linear 
combination of the functions: 

\/r J^+i/2(^V^) sin0P^(cos^) cos{ct\f\) 

V^^n+i/2(^V^) sin6> P;;(cos6>) cos(ct\/A) (36) 

where Pn is the n-degree Legendre polynomial and Jn+1/2 ^^^d 
the Bessel functions of the first and the second kind, respectively. 
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Figure 10: Two independent eigenfunctions in the case of the domain for 
A4 = As « 1.9639. 

We are concerned with the cases n — 1 and n — A. For these val- 
ues, the functions in (l36]) have, with respect to the azimuthal variable 0, 
one single bump or 4 consecutive bumps, respectively (see figure 10). As 
a matter of fact, we have: sin0P((cos0) = sin0 and sinOP'^ico^O) — 
I sin0 cos 0(7cos^ — 3). Regarding instead the radial variable r, we would 
like to find suitable linear combinations of the functions in (155]) , in order to 
impose homogeneous Dirichlet boundary conditions at r = 1 and r — td- 
In agreement with the results obtained by implementing the finite element 
method, playing with the zeros of Bessel functions, yields: 

A = A4 = A5 ^ 1.9639 and td - 5.839 (37) 

These data are now more accurate. The correct assemblage of the radial 
basis functions can be seen in figure 11. The first picture shows a linear 
combination of J3/2 and Y3/2 (^ = 1), the second one a linear combination 
of J9/2 and Y9/2 (n = 4). Up to multiplicative constants, this setting is 
unique if we impose the boundary conditions at the endpoints of the same 
interval [l,r/:)]. In a similar way, we can also update the values of table 4 
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(fourth column), namely: Ai ^ .626, A2 ^ .980, A3 ^ 1.432. 




Figure 11: Plots of combinations of Bessel functions of first and second kind, 
vanishing at r = 1 and r = td- 



Further external shells are obtained by linearly amplifying the domain D. 
In this way, the size of the successive nested shells grows geometrically by a 
factor 5.839, and the corresponding frequencies are progressively reduced by 
the same factor. Like gears of increasing magnitude connected together, the 
shells transmit their signals far away from the source, through a quantized 
process, causing a decay of the frequency at each step (see Ref. [Ij, section 
6.1). 

As a last example, we examine the case of figure 12, where the curved 
part of the domain E is an ellipse. We recall that when the Hill's vortex 
is perfectly spherical (domain A in figure 6) there is no chance to get coin- 
cident eigenvalues (see the first column of table 4). However, by reducing 
the vertical axis, the situation improves without the help of internal holes. 
We found out that A4 = A5 when the ratio between the vertical and the 
horizontal axis is equal to .6035 (see table 5). The evolution of the solution 
is expected to be similar to that of figures 8 and 9. 

In Ref. [16j, predictions and experimental verifications about the for- 
mation of stable configurations consisting of an elliptic-shaped Hill's type 
vortex, generated by an internal thin vortex ring, are discussed. Further 
indications might come from our approach, based on the study of A4 and A5, 
by simultaneously arranging the eccentricity of the ellipse and the location 
and the size of the internal vortex. The configuration is going to be similar 
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Figure 12: Elliptic section with semi-axis equal to 1 and .6035. 





domain E 


Ai 


26.96 


A2 


56.97 


As 


64.53 


A4 


101.95 


As 


101.99 



Table 5: Eigenvalues corresponding to the domain E of figure 12. 

to the one of figure 4 in Ref . [16] . We tried a few qualitative experiments in 
this direction, but, having too many degrees of freedom, a careful analysis is 
a technical exercise that we would prefer to avoid at the moment. Therefore, 
we do not anticipate any result. 

6 Conclusions 

With the help of numerical simulations, we demonstrated the possibility of 
building electromagnetic waves trapped in bounded 3-D regions of space. 
We did not insist too much on issues related to the performances of the 
algorithms. We are conscious of the fact that the methods used can be 
certainly ameliorated, in terms of costs versus accuracy. Nevertheless, this 
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was not our primary concern. 

We think that the results we got may have a general validity indepen- 
dently on the field of applications. In fact, the approach we followed, based 
on the determination of periodic solutions to the wave equation through the 
analysis of certain eigenvalues of a suitable differential operator, may be 
applied in several circumstances. For example, as an alternative, this tech- 
nique may be employed in fluid dynamics (for stable periodic flows), where 
computations are usually carried out by discretizing the equations by some 
time-advancing procedure. However, the meaning of the results obtained 
here is deeper, since they are related to the approximation of a complete 
system of hyperbolic equations (namely the Maxwell's equations) and not 
just to the detection of the flow- field. As a matter of fact, the information 
carried by our waves is not a scalar density field, but it consists of two sep- 
arated vector fields: the electric and the magnetic ones. These can be fully 
expressed by the relations in (f28ll . 

A recent subject of research, that could benefit from our investigation, 
is the detection of quantum vortex rings in superfluid helium (the literature 
in the field is very rich, see for instance Ref. [17j for a general overview). 
Note that, in liquid helium, the thickness of these rings is on the order of 
a few Angstroms. Other related topics might be the study of ball lightning 
phenomena (see for instance Ref. \18\ and Ref. [19j) or ring nebulae and 
their halos (see for instance Ref. [20]). 

Our analysis might inspire interesting applications, that at the moment 
we are unable to predict. Indeed, we showed that it is possible to detect 3-D 
regions of space, where an electromagnetic wave, suitably supplied through 
the boundary, can freely travel without bouncing off the walls. Moreover, we 
also explained how to build these regions. We are not enough experienced 
in the field of applications to find an employment for such resonant boxes, 
but this paper may suggest to some interested reader how to use them to 
create new tools or improve existing ones. 

Finally, we recall that the stability of these structures is a consequence 
of gravitational modifications of the space-time generated by the movement 
of the wave itself, according to Einstein's equation. In fact, the wave trav- 
els like a fluid, along closed stream-lines, that turn out to be geodesies of 
such a modified geometry (see also Ref. [21j). Thus, in the electrodynamic 
framework, a complete analysis of stability depends on the resolution of a 
complex relativistic problem. Nevertheless, we hope that the simple inves- 
tigation carried out in this paper may represent a step ahead towards a 
better comprehension of the mechanism of vortex formation, whatever the 
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constituents are (fields or real matter). 
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